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ABSTRACT 



^ I We present new Fokker-Planck models of the evolution of globular clusters, including 

gravitational tidal shocks. We extend our calculations beyond the core collapse by adopting 
, three-body binary heating. Effects of the shocks are included by adding the tidal shock diffusion 

' coefficients to the ordinary Fokker-Planck equation: the first order heating term, (AE), and the 

^■■^ I QApnnrl nrrlpr pnpro-^r HiQrtprQinn fprTn / A F?^ 



second order energy dispersion term, {AE ). As an example, we investigate the evolution of 



Q>«^ , models for the globular cluster NGC 6254. Using the Hipparcos proper motions, we arc now 



able to construct orbits of this cluster in the Galaxy. Tidal shocks accelerate significantly both 
core collapse and the evaporation of the cluster and shorten the destruction time from 24 Gyr 
Q I to 18 Gyr. We examine various types of adiabatic corrections and find that they are critical for 

i—^ ' accurate calculation of the evolution. Without adiabatic corrections, the destruction time of the 

C/3 i cluster is twice as short. 

^ ' We examine cluster evolution for a wide range of the concentration and tidal shock 

parameters, and determine the region of the parameter space where tidal shocks dominate the 
, evolution. We present fitting formulae for the core collapse time and the destruction time, 

H ' covering all reasonable initial conditions. In the limit of strong shocks, the typical value of the 

core collapse time decreases from 10 tj-h to 3 trh or less, while the destruction time is just twice 
that number. The effects of tidal shocks are rapidly self-limiting: as clusters lose mass and 
become more compact, the importance of the shocks diminishes. This implies that tidal shocks 
were more important in the past. 



Subject headings: stellar dynamics — globular clusters: general — globular clusters: individual 
(NGC 6254) 



Introduction 



Evolution of star clusters is driven by a variety of dynamical processes. Two-body relaxation, stellar 
evolution, evaporation of stars, and tidal truncation have all been shown to play a role. In addition to the 
internal effects, the external tidal shocks from the Galaxy are very important. Gnedin fc Ostriker (1997) 
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found that tidal shocks contribute at least as much as two-body relaxation to the destruction of the current 



Galactic sample of globular clusters. A^-body simulations of Gnedin fc Ostrikcr (1999) provide a detailed 
account of the effect of individual tidal shocks on the stellar energy distribution. In this paper we describe 
the implementation of tidal shocks in a Fokker-Planck code and present new models of globular cluster 
evolution. 



The study of the dynamical evolution of globular clusters has a long and rich history. Bpitzer (1987), 



Meylan fc Heggie (1997)| , and [Ashman fc Zcpf (1998)| provide comprehensive reviews of the evolution of 



isolated and tidally- limited clusters. Two-body relaxation causes loosely-bound stars gain velocity higher 



than the escape velocity, which leads to their evaporation from the cluster ( Ambartsumian 1938; Bpitzer 
1940). Even without evaporation, clusters would experience a catastrophic collapse driven by the negative 
heat capacity and the conductive transport of the thermal energy from the inner to the outer parts of 



the cluster (Antonov 1962; Lynden-Bell fc Wood 1968; Goodman 1987). The run-away process leads to 



a gravothermal catastrophe, or core collapse. The collapse can be reversed when a sufficient amount of 
heating is provided in the core via formation of binaries. The successive mergers of normal stars and 
subsequent supernova explosions were also shown to be effective in halting the collapse (Lee 1987). After 
the reversal of core collapse, the evolution of the cluster can be characterized by quasi-static expansion with 
possible gravothermal oscillations. 

The tidal field of the Galaxy affects cluster dynamics in two ways. The stellar distribution is effectively 



cut off at the point where the ambient Galactic density exceeds that of stars in the cluster. Lee & Ostrikcr 



(1987) found that clusters with short initial relaxation time dissolve in '-^ N/200 orbital periods around the 



Galaxy, where N is the initial number of stars in the cluster. Tidal effects may also cause cluster rotation 
as the stars on prograde orbits have higher chances to be ejected than those on retrograde orbits (e.g., 



Oh fc Lin 1992). In addition, the time- varying tidal forces cause gravitational shocks when the cluster 



passes through the disk of the Galaxy ( Ostrikcr, Spitzer, fc Chevalier 1972) or comes near the Galactic 
nucleus (e.g., Spitzer 1987). Tidal shocks increase random motion of stars and reduce cluster binding 



energy. Finally, Kundic fc Ostrikei] (1994) and Gnedin fc Ostrikcr (1997) emphasized that the tidal shock 



relaxation can be very important in accelerating cluster evolution. Recent studies (Murali fc Weinberg 
1997a-c; pstriker fc Gnedin 1997 ) show that these effects can lead to substantial evolution of the globular 



cluster population. 

All of the above processes are included in our Fokker-Planck code. We discuss the detailed 
implementation of tidal shocks and demonstrate their effects on the evolution of clusters with a wide range 
of parameters. We consider both single-mass models and models including a spectrum of stellar masses. 
Mass segregation and rapid energy transfer between the stars of different masses lead to an even faster core 
collapse and destruction of globular clusters. Although we do not include the effects of stellar evolution, 



they are likely to be important only at the early stages of cluster evolution. [Chcrnoff fc Shapiro (1987) 



and Chernoff fc Weinberg (1990) showed that the mass loss from the short-lived massive stars can disrupt 
weakly concentrated clusters (with c ^ 0.6). We also ignore the effects of rotation and velocity anisotropy 
but discuss the latter at the end of the paper. 

We first review the theory of tidal shocks in §2. The new features of the Fokker-Planck code are 
described in §3. We start with an illustrative case that demonstrates the importance of tidal shocks and 
adiabatic corrections in §4. More systematic study for a wide range of the parameter space is presented in 
§5. Finally, we discuss the implications of our results to the evolution of the Galactic star clusters. 
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2. Review of the Theory of Tidal Shocks 



The first and second order theory of the compressive (disk) shocks and tidal (bulge) shocks has been 
developed by [Ostriker, Spitzer, fc Chevaliei] (1972); ^pitzeij (1987); and [Kundic fc Ostrikei| (1994) and 



summarized in 



Gnedin fc Ostrikei (1997). Below we review the important results for the disk shocking. The 



theory of tidal shocking is discussed in detail in [Gnedin, Hcrnquist, fc Ostrikei (1999). For an alternative 



linear theory calculation see Weinberg (1994) 



In the disk shocking, the phase-averaged first and second order energy changes of stars with the initial 
energy per unit mass E are 



2 5,^ 



3^2 



Ai{x), 



E — 



9^2 



^2(2;), 



(1) 
(2) 



where is the maximum vertical gravitational acceleration produced by the disk, V is the vertical 
component of the cluster velocity with respect to the disk, and r and v are the rms position and velocity 
of stars of energy E. Here Xr,v is the position-velocity correlation factor, which takes values from —0.25 to 



0.57 (see Gnedin fc Ostriker 1999) 



The first order energy change, (AE), causes reduction in the binding energy of the system and leads to 
evaporation of the marginally bound stars. The second order change, {AE'^), causes a much larger energy 
dispersion which allows additional stars to leave the cluster. The two effects cooperate and lead to faster 
dissolution of the cluster. 

The adiabatic corrections, Ai{x) and A2{x), account for conservation of the adiabatic invariants of 
stars for which the orbital period in the cluster is short compared to the effective duration of the shock, 



H 



(3) 



where H is the characteristic scale-height of the disk. The adiabatic corrections can be approximately 
described as functions of the only dimensionless parameter, 

X = U!(r) T. (4) 

Here tu is the stellar orbital frequency of stars of energy E at the rms position r. 

Several approximations have been used to study the adiabatic corrections. (1) The impulse 
approximation is valid for the stars whose motion in the cluster is much slower than the shock, x ^ 1. This 
condition applies in the outer regions of the cluster. There both corrections become unity: 



Ai{x) ^ A2ix) ^ 1, (x<l). 



(5) 



(2) The harmonic approximation is valid for the stars which oscillate very quickly in a simple harmonic 
potential and which are not in resonance with the perturbing force. Building on results of [Spitzcr (1987), 
we get 



^is(a;) 



^2s(a:^) 



{x ^ 1; no resonances). 



(6) 



4 -t- 5e2^' 

Here we have modified the second correction from its version in Gnedin fc Ostrikei] (1997) to make it equal 
unity for a; = 0. The numeric coefficients are dictated by the asymptotic form for large values of x. (3) 
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The linear theory by Weinberg (1994)| includes the possible resonances and predicts a less steep power-law 
function. In the limit of large r, the results can be approximated ( Gnedin fc Ostriker 1999) as:|^ 



Aiw(x) = A2w(x) = (l + a;2)-3/2. 



(7) 



(4) [Gnedin fc Ostriker (1999) performed N-body simulations of the tidal shock, allowing for the self-consistent 
oscillations of the cluster potential as it relaxes into a new virial equilibrium. The simple fits to the results 
are 

^in(x) = (1 + x2)-^/2, A2nix) ^ {I + x^)-' (8) 



for the "fast shocks" whose durations are comparable to or smaller than the dynamical time at the half-mass 
radius: t <^ tdyn- For the "slow shocks", the results agree with the predictions of the linear theory (eq. [Q). 

In the case of bulge shocking, the Spitzer adiabatic corrections involve a combination of Bessel 



functions, as described in [Gnedin fc Ostriker] (1997). The Weinberg and A^-body adiabatic corrections are 
taken to be the same as for disk shocking. Since the #-body simulations of Gnedin fc Ostriker[ (1999) 
include only disk shocks, further study is needed to determine accurately the adiabatic corrections for bulge 
shocking. 

The first and second order energy changes are of comparable importance. This becomes clear if we 
define the characteristic shock timescales: 



t 



sh 



dEh/dt 

El 
dEl/dt 



\Eh\ 



disk 



disk 



{^E)h 

El 
{AE^)h 



(9a) 



(9b) 



Here Pdisk is the period of cluster's passage through the Galactic disk. The energy changes {AE)ii and 
{AE'^)h are evaluated at the half-mass radius Rh, and the characteristic energy Eh is given by the Virial 
Theorem: \Eh\ = vl^j2 « 0.2GM/Rh (e.g., [Spitzer[ 1987), where M is the mass of the cluster. In the 
impulse approximation we obtain 



l-sh 



disk 



9 p V'col 



16 



disk 



(10) 



where cuh 



s/ Rh is the rms angular velocity of stars at the half-mass radius. The two timescales are 



simply related by tsh,2 — \ tsh^ so that both processes contribute similarly to the destruction of the cluster. 



The Fokker-Planck Code 



We model the evolution of globular clusters using an orbit-averaged isotropic Fokker-Planck code 
descended from Cohn (1979, 1980). The code has been modified by Lee fc Ostriker (1987) and Lee, 



Fahlman, fc Richer (1991) to include tidal boundary and the three-body binary heating. Stars beyond the 



tidal boundary do not escape instantaneously but follow instead a continuous distribution function f{E), 
as described in Lee fc Ostriker (1987)[ This accounts for the balance of the internal and external forces at 



^This does not imply that the results of Weinberg (1994) for fast shocks are incorrect. The linear theory applies for all values 
of T, but the final expressions are complex and cannot be fitted simply. 
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the tidal radius. Therefore, those stars only drift slowly away from the cluster. The tidal field is assumed 



to be spherically symmetric, which is a weakness of the one-dimensional code (for a discussion see Lee & 



Goodman 1995). We also assume that the cluster fills its Roche lobe and require that the average density 
within the tidal radius remain constant throughout the cluster evolution. 

The heating of stars which reverses core collapse is provided by three-body binaries. They are included 



explicitly without following their actual formation and evolution, according to the prescription by Cohn 



(1985) . Ostriker (1985) argued that tidally-captured binaries are probably more dynamically important 
for massive clusters. Whether tidal capture leads to a hard binary which contributes to stellar ejection 
after a close encounter or whether it leads (more often) to a merger which causes mass loss following 
stellar evolution, the outcome is the same. Lee (1987) showed that the merger of stars would give a nearly 
identical dynamical effect on the cluster because the massive stars formed in the merger would evolve off 
rapidly. Both processes cause indirect heating by ejecting mass. The situation becomes more complicated 
if we allow for the existence of primordial binaries and massive degenerate stars. The existence of massive 
remnant stars, such as neutron stars, could cause three-body binaries to be more important than tidal 



capture binaries (Kim, Lee & Goodman 199? ). However, many aspects of the evolution after core collapse 



are independent of the actual energy source ( 



Henon 1961; Goodman 1993). 



We include the effects of tidal shocks by modifying the diffusion coefficients in the Fokker-Planck 
equation. We assume that the first and second order energy changes are known as functions of the energy 



and position using the results from 92^ and Gnedin, Hernquist, & Ostriker (1999). We now re-derive the 



Fokker-Planck equation in order to define the diffusion coefficients corresponding to the shocks. 

Let ^(E', AE) dAE be the probability of scattering of a star of energy E by the amount [AE, AE+dAE] . 



iV-body simulations ( Gnedin fc Ostriker 1999) show that the probability distribution is nearly Gaussian: 

1 



*(£;, A£;) dAE 



1 (&E-(&E))^ 

2 (AE^> 



dAE, 



(11) 



VMaW) 

where both (AE) and {AE'^) arc functions of energy E. Let N{E)dE be the number of stars in the energy 
range [E,E -(- dE]. The time evolution of N{E), in the absence of two-body relaxation and binary heating, 
can be described by the following equation 



N{E,t + At)^ J N{E - AE,t)'^{E - AE.AE) dAE, 



(12) 



where the energy changes by the amount AE in the time interval At. Expanding this equation in a series 
over AE and At, we obtain the usual form of the Fokker-Planck equation 



dt 



2dE^ 



(13) 



where {'Dt{AE))y and (23t(Ai?^))^ are the orbit-averaged diffusion coefficients of the first and second 
order. They are defined by 

_ {VtiAE))vr^dr 



mAE)), 



"""^ vr'^dr 



and 



{Vt{AE)^)y=^-° \ ^/ 



where {'Dt{AE)) and (T>t{AE^)'j are the average rates of change of the energy and its dispersion, 
respectively, per unit volume at a given location r, and r^ax is the maximum radius allowed for a star 



(14) 
(15) 
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of energy E. Given equations (P and 
function, we define 

{V,{AE'))^, ^ 



for the energy changes and equation (|llD for their distribution 



AE 
~Ai 

[AEf 
At 



^{E,AE) dAE 



{AEl 
At 



'^{E,AE) dAE = 



(AE^) + {AEY 
At 



(16) 
(17) 



The bulge and disk tidal shocks are applied at the time step when the cluster is at its perigalacticon 
or crosses the Galactic disk, respectively. The diffusion coefficients are normalized to the current time step, 
At, so as to produce the right amount of heating per event. The shocks are repeated every orbital period of 
the cluster. 

The time step was chosen such that both the cluster mass and the central density change by no more 
than 1% between time steps: 

At = 0.01 X min[(dlnp/dt)"\ [dlnM/dty^ (18) 



Lee fc Ostriker (1993) suggested that such a time step should be sufficient for an accurate calculation of 



the cluster evolution, both before and after core collapse. We compare our results without tidal shocking 



with those of Quinlan (1996), who used the same time step for the potential re-computation but broke 
it down into 32 sub-steps for updating the distribution function (the "Fokker-Planck steps"). We found 
a measurable difference in the core collapse times only for very concentrated clusters, with the initial 
concentration c ^ 2.5. For these clusters an even smaller time step may be required, although our results 
never differ by more than 30% (cf. Figure 0). 

As described below, in order to incorporate correctly the effects of tidal shocks into the code, the shock 
diffusion coefficients should be applied separately from the relaxation terms, and the potential should be 



recomputed accordingly thereafter (an alternative approach has been suggested by Johnston, Hernquist & 



Weinberg 1998). It seems awkward, though still possible, to include a set of Fokker-Planck steps for every 



potential time step in this scheme, as Quinlan (1996) did for the isolated clusters. We chose not to do so 



and therefore we had the potential updated every time the distribution function was changed. Since our 
time step (eq. JT8|) gives very good accuracy for most of the computed cluster models, we did not change it 
specifically for high concentration clusters. The resulting error should in any case be small. 



The alternative implementation of tidal shocks in a Fokker-Planck code has been done by Murali & 



Weinberg (1997a-c). Those authors employ the linear theory formalism of Weinberg (1994) and calculate 



the change of the distribution function directly from the Boltzmann equation. That approach is more 
accurate but is more cumbersome. Since the linear theory expressions are supported by the restricted 
A^-body simulations (Weinberg 1994), that method and our derived adiabatic corrections should give the 



same result. The work of Murali & Weinberg has focused more on the destruction time of globular clusters 
under various assumptions and on predicting the initial cluster distribution. In this paper, in addition 
to obtaining fitting formulae for the destruction time we present also the detailed evolution of individual 
clusters and compare it with previous work on isolated clusters. 



3.1. Comparison with A'-body simulations 

I I 

Gnedin & Ostriker (1999) report self-consistent A^-body simulations of a single disk shock of various 
strengths and durations. The cluster is modeled as a low-concentration King-model (c — 0.84) with a 
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realistic number of particles, N = 10^. The characteristic disk shocking strength is low and most of the stars 
remain bound to the cluster. However, the energy distribution changes in response to the perturbation. 
The resulting distribution after the virialization phase is accurately fitted by equations (Q), and (|^). 
Although the result depends on the initial structure of the cluster, the above fits still apply for the more 



concentrated cluster with c — 1.5 (sec Gncdin & Ostrikei 1999 for more details) 



The A'^body modeling of tidal shocks is intrinsically more accurate than the Fokker-Planck formulation. 
The former takes into account the non-linear virialization following the shock and the anisotropic effects 
induced by disk shocking. A direct comparison of the resulting energy changes between the two codes is 
difficult because of the algorithmic differences. Therefore, we have taken the simple approach to reproduce 
the A^-body results in our F-P code as closely as possible. 

The general procedure for advancing the distribution function (DP) in time consists of two parts 
(Cohn 1979). The first part introduces the first and second order changes of the DP as a solution of the 
Fokker-Planck equation. The second part updates the cluster potential by solving Poisson's equation with 
the density following from the new DF. During the second step, the DF is kept constant as a function 
of the stars' adiabatic invariants and therefore it changes as a function of energy. This procedure is no 
longer correct when treating tidal shocks which last of the order of the dynamical time of the cluster. This 
timescale is short enough that adiabatic invariants may themselves change. 

We find that in order to reproduce the A'^body results, we need to modify the above procedure for 
tidal shocks. During the time step when the cluster passes through the Galactic disk or close to the 
Galactic center, we first apply the shock diffusion coefficients and update the DF. Then we fix the DF as 
the function of energy and obtain the new potential by solving Poisson's equation. Specifically, since the 
potential is changed and the energy grid is renormalized, we require that in each energy bin the new energy 
distribution, Nnew, scale with the new energy grid, Enew ^ E + AE, as 

iV„e«,(S«e») dE„,^ = N{E) dE, (19) 

where N{E) is the energy distribution before the potential recalculation. We have checked that the above 
procedure reproduces well the A^-body results for a range of shock amplitudes. Figure |l| shows the change 
in the energy distribution after a single shock in a test cluster. The agreement between the Fokker-Planck 
and A/^body results is reasonably good. Also, Figure ^ shows that the previous procedure for updating the 
DF for the fixed adiabatic invariants departs much more strongly from the A'^body results than docs our 
new prescription. 

Another effect demonstrated by the A^-body simulations is a slight decrease in the shock-induced energy 
dispersion subsequent to and caused by the potential fluctuations following the shock. This paradoxical 



"refrigeration effect" persistently survived all our tests and seems to be real (see Gnedin & Ostriker 1999). 
However, the magnitude of the dispersion decrease is relatively small, which led us to ignore it in the 
present calculation. Once the nature of this effect is understood better, it should also be included in the 
models of globular clusters. 



4. Typical Evolution of a Single Cluster 

We consider now the detailed models of the globular cluster evolution. For the most part, we study 
single-mass models and then briefly investigate the effect of a spectrum of stellar masses. 

With the new data on proper motions from Hipparcos, it is now possible to reconstruct full three- 
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dimensional velocities for some of the Galactic globular clusters (Odenkirclien et al. 1997) and to build 



realistic models of their evolution. We use the inferred orbital parameters for NGC 6254 and calculate its 



orbits in an analytic potential of the Galaxy from Allen fc Santillan (1991) . The orbits are of the rosette 



type (Figure ^ with a small eccentricity, e — 0.21. The cluster moves around the Galaxy with a period of 
about 1.4 X 10* yr, during which it passes next to the Galactic center every 9.5 x 10^ yr and crosses the 
Galactic disk every 5.3 x 10^ yr. The period ratio differs from the expected factor of two due to the slow 
precession of the orbit. 

From the computed orbits of NGC 6254 we infer that the amplitude of tidal shocks varies approximately 
as a Gaussian function of time. This satisfies the conditions under which the A^-body adiabatic corrections 



were calculated (see Gnedin & Ostriker 1999) 



4.1. Single Mass Models 



The observed structure of NGC 6254 is well characterized by a King model (King 1966) with the 
concentration c = 1.4, or the structural parameter Wq = 6.55. The mass of the cluster is 2.25 x 10^ -^0? 
assuming a constant mass-to-light ratio M/Ly = 3 in solar units. First, we consider in detail a model 
composed of stars of the same mass, = 0.7 Mq. The important parameters of the model are summarized 
in Table |l|. 

We implement both disk and bulge shocks in our integration. The models differ in the way the shocks 
are included: no tidal shocks, steady tidal field and relaxation only {case 0); tidal shocks with the first 
order (AE) term only {case 1); tidal shocks with the second order {AE^) term only {case 2); and finally 
tidal shocks with both (AE) and {AE"^) terms {case 3). 

An important observable of the cluster evolution is the mass-loss due to tidal shocks and due to the 
direct escape of high- velocity stars. Figure ^ shows the run of the cluster mass with time. As is well-known, 
the mass of a tidally limited cluster goes to zero in an approximately linear fashion. A convenient choice is 
to express time t in the units of the initial half-mass relaxation time trh.o (eq. l2l|| ). Two-body relaxation 
leads to destruction of the cluster in about 32trh,o- Note that the numerical problems cause the cluster to 
lose its stability and dissolve into the sea of background stars before its mass reaches exactly zero. Usually 
the code fails to re-calculate the cluster potential when the mass falls to about 1% - 4% of the initial 
mass. However, we can extrapolate through the last few time steps in order to estimate the hypothetical 
time when the cluster mass vanishes. The disruption time obtained this way is an overestimate of the real 
disruption time, but all other estimates would suffer from a personal choice. Thus, in the rest of the paper 
we use the extrapolation of our calculations as the destruction time td- 

The destruction time is significantly reduced when we include gravitational shocks. The {AE) term 
alone changes td from 32tr/i.o to 2(!>trh,o and the shock-induced relaxation brings it to about 2Atrh,o- To 
access the importance of the two shock terms separately, we turn off the energy shift, (Ai?), for the case 2 
run. The second order term leads to an enhanced mass loss relative to the relaxation case for most of the 
evolution. This effect is not as strong as the first order effect due to the limiting adiabatic corrections, which 
we investigate in more detail laterQ. For comparison, we also show the mass loss that would have been due 
to the shock-induced relaxation if there were no adiabatic corrections (the impulse approximation). The 
latter effect would reduce the destruction time of the cluster by about a half. 



^The effect of the second order term is reduced relative to the results of 



Gnedin & Ostriker 



(1997) by inclusion of the 
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An important feature of tidal shocks is the enhanced mass loss at the early stages of the evolution, 
prior to core collapse. In the pure relaxation case, the mass loss is slow until core collapse and linear in time 
afterwards. The first theoretical calculations of the probability of stars' escape from the cluster through 
two-body relaxation ( Ambartsumian 1938; ^pitzci] 1940) predicted the following dimcnsionlcss rate 



trhjt) dM 
M(t) dt 



0.0074, 



where trh is the half-mass relaxation time (Spitzer & Hart 1971 



Ml/2 ^3/2 

0.138- 



61/2 ln(A) ' 



(20) 



(21) 



Here M is the current cluster mass, Rh is the half-mass radius, m* is the average stellar mass, and 
ln(A) = ln{OAN) is the Coulomb logarithm, being the number of stars in the cluster. Later calculations 
of the tidally-truncated cluster by Henon (1961) gave a larger value, ~ 0.045. 

I I 

Spitzer & Chevalier (1973) used Monte-Carlo simulations for several models of globular clusters and 

found = 0.05 for Rt/Rh = 3.1, and = 0.015 for Rt/Rh — 9.3, where Rt is the tidal radius of the cluster. 
Based upon those simulations, Aguilar, Hut, fc Ostriker (1988; hereafter AHOj proposed the following 
fitting formula for the time to disruption: 



td = 0.15 



Rh 
Rt 



t 



rh,0: 



(22) 



where trh,o is the initial relaxation time. For NGC 6254, the disruption time predicted by AHO would be 
about 52trh,o^ whereas we see that the cluster dissolves in 32ir/i,o even without tidal shocks. In general, we 
find that equation ( |2^ ) overestimates the actual disruption time by a factor of several. 

We compare the previous results for the escape probability with our Fokker-Planck calculations in 
Figure IJ. In the relaxation model {case 0), the value of rises almost monotonically through core collapse 
and until destruction of the cluster, reaching and exceeding Henon's estimate only in the late stages of 
the evolution. On the contrary, in the shock-dominated models the escape probability is very high in the 
beginning when tidal shocks efficiently remove stars from the outer parts of the cluster. The mass loss then 
slows down and conforms to the relaxation rate soon before core collapse. Then again, rises with time to 
reach Henon's self-similar value. This later increase caused primarily by the declining mass of the cluster, 
since the mass loss rate is almost constant in the post core-collapse cluster. 

The cluster structure can be described in terms of the characteristic radii, such as the core radius, 
reflecting compactness of the cluster; the tidal radius, confining all stars bound to the cluster; and the 
half-mass radius, relating to the global properties of the cluster. Figure ^ shows the evolution of these 
parameters normalized to the initial core radius. The core radii for all models drop to extremely small 
values at the point of core collapse. The core-collapse time tec is about 13 trh.o for the relaxation case, in 



agreement with Lee fc Ostriker (1987) and Quinlan (1996) . Tidal shocks speed up core collapse significantly: 
tec = 10 trh,o in case 1 and it is still smaller when we include the shock-induced relaxation term. The 
half-mass radius changes only by a factor of two over the whole cluster lifetime, which indicates that core 



correlation factor ( 1 -|- Xr.v) (eq. At each time step, we calculate the value of Xr.v{E\ c) using the fitting formula from 

Gnedin fc Ostriker (1999). 
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collapse affects only a small fraction of stars. The tidal radius decreases slowly, reflecting the gradual mass 
loss from the cluster. 

Figure ^ illustrates core collapse more clearly. The central density of the cluster, rises by nine 
orders of magnitude before the collapse is halted by a production of three-body binaries in the very dense 
core of the cluster. The energy released from the binaries eventually reverses the collapse. The cluster 
consequently expands, though the central density still remains very high (about four orders of magnitude 
higher than initially). The actual value of pc at the core bounce and post-collapse phase depends on the 
specific heating mechanism. A stronger heating rate gives a lower central density during the post-collapse 
phase. 

The cluster concentration, c = log]^o(^t/^c), is closely related to the central density. Since the tidal 
radius Rt decreases slowly with time as M^/^, the concentration varies roughly as logpc- 

The core collapse time provides a very useful yardstick for cluster evolution. Even in the presence of 
tidal shocks, the evolution of the models is similar when the units of time are scaled to the core collapse 
time. Figure illustrates this point. The central density in models 0-3 varies almost identically with t/tcc- 
Therefore, we can use this scaling to parameterize cluster evolution. The destruction time is an almost 
constant multiple of the core collapse time: td ~ 2.5 tec- 

The post-collapse evolution of globular clusters depends on the concentration and the remaining mass. 



The self-similar solution by Henon (1961) for tidally-limited clusters can be expressed (Goodman 1993) as 
follows: 

U-t ^22 Atrh{t), (23) 

where — t is the time to disruption at time t. Figure ^ shows that the self-similar limit is approached 
at late stages of the cluster evolution, but for the most part equation (^3|) underestimates the time to 
destruction by about 50%. 



4-.1.1. Comparison of adiabatic corrections 

An important constituent of our models is adiabatic corrections for the tidal shock effects (eqs. ^-||)- 
Various corrections predict different impact of shocks in the middle parts of clusters and lead to different 
evolutionary paths. The amplitude of the corrections depend on the ratio of the perturbation timescale to 
the dynamical time of stars in the cluster. For NGC 6254, the characteristic time for disk shocking (eq. j|]) 
is Tdisk ~ 1-3 X 10^ yr, and for bulge shocking (cf. [Gncdin, Hernquist, fc Ostrikei 1999) is Tbuigc = 1.3 x 10'' 



yr. The half-mass dynamical time of the cluster is initially much shorter, ^ =3.1 x 10^ yr, leading to 
strong suppression of the effects of tidal shocks. 

Figure ^ compares the Spitzer, Weinberg, and iV-body adiabatic corrections for the case 3 models of 
NGC 6254. When the shocks are "slow" (a; = iVhT ^ 1), as indicated by the true values of Tdisk and Tbuigo, 
both SI and Nl models have tidal shock effects suppressed more strongly than does the Wl model with the 
less steep adiabatic corrections. For the "slow" shocks, the Wl model gives the correct description of the 
cluster evolution. 

When we decrease the shock timescales by a factor of five, we enter the regime of "fast shocks" and 
the N2 model is correct. For a; ^ 1, we can estimate the importance of the adiabatic corrections by the 
Taylor expansion of equations (^-^. The N2 model shows the slowest evolution because Ai^ « 1 — 2.5x^, 
whereas for the W2 model, Aiw ~ 1 — 1.5a;^. Even though the first order correction for the S2 model. 
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Ais « 1 — 2x'^, is smaller than Aiw, the second order correction is larger, A2S ~ 1 — 1.1.t^. This allows the 
S2 model to evolve faster than the other two. 

Finally, when we reduce the shock timescales by a factor of hundred, the impulse approximation applies 
everywhere in the cluster and all three models are essentially identical. The evolution in this case proceeds 
much faster: the cluster is destroyed in half the time of the original model Wl. Thus, adiabatic corrections 
are critical for calculating correct globular cluster models. 



4- 1.2. Low concentration model 

For comparison, we consider another model of NGC 6254, with lower concentration c = 0.84, or 
Wq ~ 4. This is the model used for A'^-body simulations of Gncdin fc Ostriker| (1999). All other parameters 



are the same as in the previous model. We keep the observed tidal radius in parsecs, so the core radius and 
the relaxation time scale in physical units accordingly (see Table ||). Therefore, both models have the same 
average density, pav oc M/Rf. 

Figure |l^ shows the mass evolution of the low concentration model. The cluster evaporates much 
faster than the high concentration model, in units of the initial relaxation time trh,o- Clearly, tidal shocks 
have much stronger effect on less concentrated clusters but the total speed-up in time to destruction, due 
to shocks, is again about 25%. Notice that the destruction time in years is about the same for the low- and 
high-concentration models. Figure |ll| shows that the evolution of the central density is qualitatively similar 
to the previous model. However, the effect of the tidal shock relaxation is stronger for the low-concentration 
model. 



4.2. Multi-mass models 

In a single mass model, the energy is transferred between the cluster core and the envelope, which 
requires a conduction mechanism acting on the half-mass relaxation timescale. In a multi-mass model, the 
energy exchange occurs also between stars of different mass. A typical timescale for the energy exchange 
between the two mass species mi and 1712 (where 7712 < mi) is trh x (to2/toi). Since this is much shorter 
than the core collapse time for single mass clusters, the overall evolution of multi-mass clusters is expected 
to be faster. 
I 1 

Lee & Goodman (1995), for example, considered the evaporation of a multi-mass cluster in a steady 
tidal field and found that the mass loss rate can more than double relative to the single-component case. 



Chernoff & Weinberg (1990) investigated the effects of stellar evolution on the cluster dynamics. In the 
future, it will be important to generalize the present calculations by including a realistic mass function and 
stellar evolution, along with the presently implemented two-body relaxation and tidal shocks. In this paper, 
we present only an illustrative model allowing for the mass spectrum. 

We take the previously described model for NGC 6254 and include seven mass components, ranging 
from O.IMq to O.TMq. The number of stars in each component N{m) oc m~^. Thus, the upper mass limit 
corresponds to the turn-off of the main sequence at 10 billion years (we also used it in the single-mass 
models), and the lower limit is a reasonable cutoff of the initial mass function. The relaxation time is 
defined by equation (pT]), where we substitute for the mean stellar mass (see Lee fc Goodman 1995). 
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Figure 12 shows the mass loss for the multi-mass model. Notice that although the evolution proceeds 
faster in units of the initial half-mass relaxation time, this timescale expressed in years is longer than that 
for the single-component model (Table Q). As a result, the destruction time is roughly 20 Gyr in both cases. 
The contribution of tidal shocks in the multi-mass model is greater, reducing the time to destruction by 
about 37%. 



4-2.1. Comparison with N-body models 



Recent paper by Takahashi & Portegies Zwart (1998) compares the multi-mass two-dimensional 



Fokker-Planck (F-P) models of Takahashi (1995) with the A^-body simulations using a special purpose 
computer, GRAPE-4-. The models include stellar evolution, two-body relaxation, and the tidal truncation 



for a low concentration (c — 0.67) cluster model from Chernoff & Weinberg (1990). The above authors 
argue that the escape of stars through the tidal boundary in an isotropic F-P code, where the distribution 
function depends only on the energy, is faster than in the A'^body models. To remedy this disagreement, 
they suggest a new escape criterion based both on the energy and angular momentum of stars. This 



■'apocenter" criterion can only be implemented in the anisotropic F-P code of Takahashi, Lee, fc Inagaki 



(1997) 



The old isotropic and the new anisotropic escape criteria lead to significantly different results only 
when the mass loss rate per relaxation time is high. This especially affects the early stages of the cluster 
evolution, when stellar winds from young massive stars drive mass loss on a relatively short timescale. 



Immediate removal of that mass in the Chernoff & Weinberg (1990) models lead to the disagreement with 
the A^-body simulations (see Takahashi & Portegies Zwart 1998). Since our models do not include stellar 
evolution, we expect less disagreement. 



We run a case model of the cluster using the parameters from phcrnoff fc Weinberg (1990) . We 
use 14 mass components from O.AMq to 15Mq, distributed according to the specified mass function, 
N(m) cx m~^'^. Without stellar evolution the time to destruction in our model should be longer than in the 
A'^body and the anisotropic F-P models which give td ~ 1.3trh,Q. Indeed, we find that our model evolves 
slower, with « 2trh,o- This indicates that our results represent an upper limit on the destruction times of 
globular clusters. 



5. Review of Globular Cluster Evolution Including Tidal Shocks 

The evolution of the isolated single-mass King models (in dimensionless units) is determined by the 
only parameter, the concentration c. Inclusion of tidal shocks introduces another independent parameter, 
for example, the ratio of the half-mass relaxation time to the characteristic shocking time, 

P = trh/tsh- (24) 

The parameter /3 is a combination of two variables describing the structure of clusters: the number of 
stars, N, and the characteristic density, ph oc M/R\. When the amplitude of the tidal force on the cluster 
is fixed, the shocking time (eq. (l0|) scales as tsh « ph, and the relaxation time (eq. [0) scales as 
Uh oc p^^/^A^/lnA, where A = 0.4A^. Then 

r, N _3/2 , s 

P^^Pn'- (25) 
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While the concentration c describes the effects of two-body relaxation, the parameter (3 shows the relative 
importance of tidal shocks. In this Section, we explore the evolution of globular clusters for a wide range of 
initial parameters c and /3. 

All models in this section include fully the relaxation and tidal shock effects (i.e., case 3) with the 
adiabatic corrections for "slow shocks", equation (Q). We arrange a grid of initial conditions, covering 
a number of cluster concentrations, from c — 0.6 to c = 2.6 with the step 0.2, and a range of shock 
parameters, (3. We compute 15 models per concentration family, equally spaced on the logarithmic scale 
from (3 « 10~^ to /? « 10^. All other parameters correspond to the model of NGC 6254 described in the 
previous section. The observed tidal radius of NGC 6254 is fixed and the core radius scales accordingly 
with the concentration. Thus all models have initially the same average density, M/Rf. 



5.1. Cluster evolution in figures 

Figure |l^ shows how different regimes of cluster evolution map on the plane of parameters c and /3. In 
the lower left region of the plot tidal shocks are weak and evolution is dominated by two-body relaxation; 
in the upper left region core collapse proceeds on a relatively short timescale without affecting most of the 
cluster; and in the lower right region tidal shocks are strong and cause substantial mass loss before core 
collapse. The upper right region is never reached in reality because of extremely strong tidal shocks leading 
to fast destruction of the clusters. 

Qualitative evolution in the left region of the plane is characteristic of a tidally-truncated model. As 
clusters start to collapse, the relaxation time becomes increasingly shorter and the tidal shock time becomes 
increasingly longer because the mean density ph rises (Figure ^4|). The mass loss is weak at this stage and 
the number of stars, N, stays relatively constant. As a result, the parameter (3 decreases (and moves to the 
left of the diagram) until the very advanced stages of core collapse. At that point the value of (3 freezes and 
the evolutionary tracks are strictly vertical (the upper part of the diagram). After core collapse, the cluster 
re-expands and its density falls slightly, moving it to the right of the c — (3 diagram. At a later time, so 
little mass is left that two-body relaxation speeds up again and drives clusters to complete dissolution. 

Figure |l^ illustrates the behavior of the parameter /3 with time. On the low end, all lines scale 
similarly, following the evolution of the mean density. For the large initial values of f3, the early mass loss 
due to shocks is very important. This leads us to investigate the right part of the c — (3 diagram. 

The evolution of clusters with strong tidal shocks differs dramatically from the previous case. An early 
mass loss caused by the shocks changes the structure of the clusters, removing stars from the outer parts 
and adding energy dispersion in the core. This increases the central and the mean density of the clusters. 
Both core collapse and final destruction proceed much faster. These effects of tidal shocks leave noticeable 
ripples in the clusters tracks on the right part of the diagram (Figure |l3|) . 

Tidal shocking is rapidly self-limiting. Clusters with large values of f3 quickly lose mass and evolve 
to /3 ^ 0.1. Figure |l^ demonstrates that just a first few shocks lower the value of /3 by several orders 
of magnitude. Every subsequent shock causes the see-saw variations of the shock parameter, with an 
increasing amplitude towards the late stage of evolution when fewer stars are left. 

A few clusters in the right part of the c — [3 diagram cross their evolutionary tracks on the way to core 
collapse. At the time of crossing, at least one of the clusters has already suffered a severe mass loss caused 
by the shocks and has a structure significantly deviating from King models. The density profiles of the two 
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clusters are similar within the half- mass radii but depart from each other closer to the tidal radius. Those 
clusters with higher initial concentration and stronger shocks are more severely truncated than the clusters 
with smaller concentration and weaker shocks. Also, tidal shocks add the velocity dispersion in the core 
and, therefore, raise the central density so that the ratio Rt/R^, as measured by the concentration c, is the 
same for both clusters. 

At the crossing point the two parameters, c and /3, no longer uniquely specify cluster's future path. A 
third parameter comes into play as we investigate the thermodynamics of clusters. 



5.2. Cluster thermodynamics 

From a thermodynamic point of view, we can consider stars as particles of gas. A critical factor 
determining evolution is the heat flow between the core and the halo of the cluster. The heat in this case 



is the velocity dispersion of stars, T = w^, whose flux is derived in |Spitzer| (1987, p. 68). The amount of 
heat transported through the cluster per unit time, or the "luminosity", is 47rr^_F'(r). The heat conduction 
facilitates core collapse and proceeds on a relaxation timescale. Let us construct a new dimensionless ratio 
of the heat transported through the half-mass radius in the relaxation time trh, to the total reservoir of 
heat, T{rh)- Neglecting constant factors of order unity, we define the following parameter 

d In 



1 



dlnr 



(26) 



which is essentially a logarithmic gradient of the velocity dispersion. For an isothermal sphere, 7 = 0. 
In general, for clusters on the way to core collapse, this parameter is positive as the heat is transferred 
from the kinematically hot core to the cold halo and grows as core collapse accelerates. In contrast, tidal 
shocks try to reverse the heat flow by heating preferentially the outer parts of the cluster. Even though the 
instantaneous value of 7 is determined only by the density structure, its derivative, 7, is a signature of the 
relaxation- or the tidal shock-driven evolution. 

The clusters crossing paths on the c — (3 diagram have similar values of 7, but the derivatives have 
opposite sign. For the clusters moving almost vertically on the diagram 7 > 0, whereas for the clusters 
which cross the diagram horizontally 7 < 0. Therefore, the evolution of the former is dominated by 
two-body relaxation and of the latter by tidal shocks. 

Knowing the value of 7 for the Galactic globular clusters, we could in principle differentiate the two 
evolutionary paths. Unfortunately, it is hard to establish the sign of 7 observationally. 

Another signature of the tidal shock-dominated clusters is a steep density profile and a low number 
density of stars at the tidal radius. Strong shocks sweep away a large number of stars at once and leave the 
cluster filling smaller volume in space than the Roche lobe imposed by the external tidal field. Clusters 
with a sharp density contrast and a high velocity dispersion of the halo stars might be experiencing strong 
shocking. Most of such clusters are expected to come close to the Galactic center but are not necessarily 
present there now. 
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5.3. Cluster evolution in numbers 

The mass loss from clusters is strongly enhanced by tidal shocks. An illustration of this is the fraction 
of the initial mass remaining at the time of core collapse, M{tcc)- Figure |l^ shows that the larger /?, the 
smaller the remaining mass. Note, that many models of various initial concentration converge at /5 « 0.1, 
giving M(tcc) ~ 45%. 

The core collapse time is a strong function of both initial parameters c and /3. Figure ^ shows tec for 
the models without tidal shocks. In units of the initial half-mass relaxation time, the core collapse time can 
be fitted as follows: 

^'^^ _ /l(c) — \Q''-^+'^^'^+''-3C^+a4C^+a5c'^ (27) 
trhM 

where the coefficients ai are given in Table |^. Our results are in good agreement with iQuinlan (19961 



except for the clusters with very high concentration. Even for them, the maximum difference in tec is less 
than 30%. 

Figure |l^ shows the core collapse time for the models including tidal shocks. The shock-induced 
relaxation speeds up core collapse by adding velocity dispersion in the core and by removing stars in the 
outer parts of the clusters, thus reducing the relaxation time. We fit the results assuming, as an Ansatz., the 
following functional form: 

trh.O trh,0 



~tj Cf if 



(1 + 61/3^^+63/?''"). (28) 

13=0 

The coefficients bi for models of various concentration c are given in Table |^. The fit is generally accurate 
to 20%, with larger errors for models with very high initial concentration. 

The time to destruction, td, is also strongly affected by tidal shocks. As we have emphasized in the 
previous section, evolution of clusters with different parameters scales similarly when time is normalized to 
the core collapse time. Accordingly, we provide the fitting formulae for the scaled destruction time, td/tcc- 



For the models without tidal shocks (Figure 19), we use the same functional form as equation (E7|) and fit 



td both in units of the initial relaxation time and in units of the core collapse time: 



/2(c), :?^ = /3(c), (29) 



trh,0 t 

where the coefficients ai are again given in Table || 

For the models including tidal shocks (Figure the scaled destruction time varies very little 
for the low concentration clusters and increasingly more for the high concentration ones. Interestingly 
enough, regardless of the initial structure, models with strong shocks converge to an almost constant value 
td/tcc ^ 2. It follows that the core collapse time is a midpoint in cluster's lifetime. 

Since the value of td/tcc is almost constant for /3 1 and for /3 ^ 1, we choose the following functional 
form to fit the results: 



td _ td 

tec tec 



1 + dif]'^^ 



\ ^^rl R " (30) 

where the coefficients ai are given in Table IJ. This expression allows for a variation in the slope of the fit 
for large values of /3, but the coefficient c?2 is always close to unity, as expected. 

In the case of very strong shocks, both the core collapse and the destruction times scale approximately 
as tee, td oc I3~^/'^ (cf. the coefficient 62 in Table 0). This asymptotic behavior can be understood as follows. 
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The shock parameter /3 is proportional to the mean energy change due to shocks (eq. which in turn is 
proportional to the amplitude of the tidal force squared, . Therefore, the approximate scaling of tec and 
td indicate that these timescales vary inversely proportional to /. 



6. Discussion 

Our Fokker-Planck models provide the most comprehensive calculations of the globular cluster 
evolution including gravitational tidal shocks. We find that tidal shocks significantly alter the evolution. 



in agreement with the earlier studies by Spitzer fc Chevalier (1973) and more recent ones by Gnedin & 



□striker (1997) and Murali & Weinberg (1997a-c). Globular clusters appear to be rather fragile with regard 



to the effects of tidal shocking. 

Tidal shocks accelerate dynamical evolution and evaporation of clusters by effectively stripping stars 
in the outer regions and adding energy dispersion in the core. These effects are self-limiting: as clusters 
lose mass and become more compact, the relative importance of tidal shocks rapidly diminishes in favor 
of two-body relaxation. In the case of the Galactic globular cluster NGC 6254, tidal shocks decrease the 
destruction time from 24 Gyr to 18 Gyr. 

We examine several types of adiabatic corrections: the harmonic approximation (of L. Spitzer), the 
linear perturbation theory (of M. Weinberg), and the results of the A^-body simulations. The differences 
between the various types of adiabatic corrections depend on the regime of "slow shocks" or "fast shocks" , 
but overall the corrections are critical for accurate calculation of the evolution. Without adiabatic 
corrections, the destruction time of NGC 6254 is twice as short. 

The evolution of models with a single mass component can be characterized by the King-model 
concentration parameter, c, and the shock parameter, [3 (eq. [Q). The effects of tidal shocks are important 
when j3 ^ 10~^ — 10~^, depending on the concentration. Tidal shocks increase both the mean and the 
central densities, relative to the relaxation case. Core collapse is accelerated by about 30% for (3 — 0.01 
and by a factor of 3 for /3 = 1. In the limit of strong shocks, the core collapse time scales approximately 
as tec oc f3~^/'^. The time to destruction is also dramatically reduced, following the same scaling law. 
Overall, we find that evolution of models proceeds very similarly when time is expressed in units of the core 
collapse time. In the limit of strong shocks, the destruction time td ~ 2tcc, independently on the initial 
concentration. Core collapse provides a useful midpoint of cluster evolution. 

We also consider an illustrative model with a range of stellar masses. In units of the initial relaxation 
time, the evolution proceeds much faster than in the single-mass case, but for both models the time to 
destruction is comparable when expressed in years. Thus, the models with the same average density evolve 
in approximately the same time. A systematic study of multi-mass models is beyond the scope of the 
present study and requires a knowledge of the initial mass function in star clusters. 

The results of our calculations have important implications for the evolution of globular clusters. For 



the sample of the Galactic globular clusters from Gnedin & Ostriker (1997) the typical values of the shock 
parameter are /3 ^ 10"^ — 10^^. This indicates that the effects of tidal shocks and of two-body relaxation 
are comparable at present. Since the shock effects quickly saturate, tidal shocks must have been more 
important in the past. Not only will the current population of clusters be significantly depleted in the 



future ( Gnedin fc Ostriker 1997), but also a large fraction of the initial population might have been already 



destroyed in the Galaxy (Murah & Weinberg 1997c) 
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Present calculations can be improved in several ways. A reasonable mass function and effects of stellar 
evolution should be included in order to accurately describe the early phase of cluster evolution. Also, in our 
treatment of tidal shocks we have neglected small (but real) negative relaxation due to cluster oscillations. 
The final weakness of our models is the assumption of the isotropic velocity distribution. Radial anisotropy 
would develop in the outer parts in the absence of tidal shocks. Since the stars on radial orbits will be 
ejected more easily than those on circular orbits, tidal shocks should effectively isotropize orbits even in 
the outer parts. Takahashj (1995,1997) and Takahashi, Lee, & Inagaki (1997) have developed an accurate 



method for integrating the anisotropic Fokkcr-Planck equation. We plan to include the effect of tidal shocks 
in the anisotropic code in order to clarify these points. 

We would like to thank Kathryn Johnston for useful discussions, and the anonymous referee and the 
scientific editor Steven Shore for detailed comments. This work was supported in part by NSF grant AST 
94-24416 and by Korea Science and Engineering Foundation grant No. 95-0702-01-01-3. 
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Table 1. Model parameters for NGC 6254. 



Model c trhfi (Gyr) tsh/Uhfi tcc/trh,o td/Uhfi 



Single mass 1.4 0.76 96 

case 12.9 32 

case 1 10.0 26 

case 2 11.2 32 

case 3 9.7 24 

Low CONCENTRATION 0.84 1.80 7.3 

case 12.2 16.5 

case 1 8.0 13.5 

case 2 8.5 14.5 

case 3 7.1 12 

Multi-mass 1.4 2.77 15 

case 1.76 13.5 

case 1 1.59 9.5 

case 2 1.60 13 

case 3 1.57 8.5 



Table 2. Fitting coefficients for tec and td versus c. 



Function ai a2 as 04 05 



/i(c) 1.3162 -1.9809 3.2910 -1.8401 0.2987 

/2(c) 1.2897 -1.1223 1.9824 -0.8803 0.1104 

/3(c) 0.1286 0.3173 -0.6573 0.6367 -0.1319 
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Table 3. Fitting coefBcients for tec versus /3. 



c 


h 


b2 


b3 


b4 


0.8 


1.5029 


0.5539 


0.0185 


0.4558 


1.0 


0.1174 


1.4184 


1.6004 


0.4491 


1.2 


0.5124 


0.8597 


1.5189 


0.3840 


1.4 


3.3268 


0.4680 


-1.0819 


0.5281 


1.6 


2.9924 


0.4613 


-1.0366 


0.5726 


1.8 


1.1447 


0.3691 


0.1791 


0.3614 


2.0 


1.4912 


0.2996 


-0.7153 


0.1606 


2.2 


0.8931 


0.6518 


-0.6721 


0.2626 


2.4 


0.6413 


0.5716 


-0.9622 


0.2167 


2.6 


0.9452 


0.5916 


-1.4586 


0.2065 



Table 4. Fitting coefficients for ta versus (3. 



c 


di 


d2 


d3 


0.8 


2.1942 


1.0376 


2.3645 


1.0 


0.1315 


1.2982 


0.2646 


1.2 


1.6051 


1.1528 


2.6816 


1.4 


45.420 


1.0356 


74.731 


1.6 


135.52 


0.9993 


300.84 


1.8 


188.62 


0.9932 


587.29 


2.0 


159.24 


0.9675 
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Fig. 1. — Comparison of the change of the energy density, AN{E), for the current Fokker-Planck models 
(dashes) and the A/^body calculations [solid lines). Energy on the horizontal axis is normalized to the cluster 
binding energy. The test cluster is initially a King model with the concentration c = 0.84. A single weak 
shock is applied, similar to the ones studied by Gnedin & Ostrikei (1999). Top panel: the energy change 



due to the shock; Bottom panel: the energy change due to the potential readjustment. Dots show previous 
implementation of tidal shocks in the F-P code, where the distribution function was kept fixed as a function 
of adiabatic invariants. The current procedure agrees with the A'^body results more accurately, especially in 
the center and in the outer parts of the cluster. 
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Fig. 2. — Orbits of NGC 6254 in the analytic potential of the Galaxy from 
period of 3 x 10^ years. The dot marks the current position of the cluster. 
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Fig. 3. — The mass loss for the single-mass models of globular cluster NGC 6254. Time is expressed in units 
of the initial half-mass relaxation time, as well as in billions of years. Dots are for the model with two-body 
relaxation only {case 0), dashes are for the model that includes the effects of the energy shift due to tidal 
shocks (case 1), long dashes are for the model that includes the second order energy dispersion (without the 
heating term; case 2), and the solid line is for the final model including proper treatment of all of the shock 
effects {case 3). The adiabatic corrections for "slow shocks" are used (eq. 0). For illustration, thin dashes 
show the effect of the second order term (marked "^imp") if it had not had adiabatic corrections. Filled 
circles indicate the time of maximum core collapse. 
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Fig. 4. — Escape probability in units of the half-mass relaxation time for NGC 6254. Line notation is 
the same as in Fig. ^j. The lower horizontal line shows the analytical estimate by Ambartsumian (1938), 
= 0.0074, and the upper self-similar calculations by Hcnorj (1961), — 0.045. 




Fig. 5. — Evolution of the tidal, half-mass, and core radii of NGC 6254. Line notation is the same as in Fig. 
i 




Fig. 6. — Evolution of the central density of NGC 6254. Line notation is the same as in Fig. ^ 
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Fig. 8. — Time to destruction of the cluster NGC 6254 in units of the current relaxation time, trh{t), vs. 
cluster age. Line notation is the same as in Fig. ^ The filled circles indicate the time of maximum core 
collapse. The dotted horizontal line is Henon's self-similar solution, equation (^). The abscissa is normalized 
to the different destruction time td for each model. 
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Fig. 9. — Comparison of the Spitzer {S), Weinberg ( W), and A^-body {N) adiabatic corrections for the case 3 
models of NGC 6254. The first set of models {SI, W1,N1) is calculated with the true values of the disk and 
bulge shock timescales, Tdisk and Tbuigc; the second set {S2,W2,N2) with the shock timescales reduced by a 
factor of 5; and the third set {S3, W3,N3) with the shock timescales reduced by a factor of 100. Solid lines 
correspond to the Weinberg corrections (eq. 0), dots to the Spitzer corrections (eq. [^), and dashes to the 
A^-body corrections for fast shocks (eq. [^). Filled circles indicate the time of maximum core collapse. 




Fig. 10. — The mass loss of the low concentration model, c=0.84. All other parameters of NGC 6254 are 
fixed, including the tidal radius in parsecs. Line notation is the same as in Fig. ^. 




Fig. 12. — The mass loss of the muhi-mass model of NGC 6254. Line notation is the same as in Fig. ^. 
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Fig. 13. — Evolution of the model elusters as a function of the concentration, c, and the shock parameter, 
/3. Stars show the initial conditions, filled circles the time of core collapse. Marked regions on the diagram 
correspond to the different regimes of cluster evolution. The critical values of /?, separating the left and right 
regions, are determined when the core collapse time deviates from its value for /3 = 0. 
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Fig. 14. — The mean density, ph, as a function of time normalized to the individual core collapse time for 
each run. All models shown have the same initial concentration c = 1.4 but different shock parameters 
varying from (3 = 10~^ to /3 = 10. Models with (3 < 10~^ show essentially identical evolution. At the time 
of core collapse, the mean density is higher for the higher values of f3. 
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Fig. 15. — The shock parameter, /?, as a function of time normahzed to the individual core collapse time, 
for the models with the various initial values of (3 and the same initial concentration c = 1.4. 




Fig. 16. — The mass remaining at the time of core collapse vs. the shock parameter, /3, for the various 
families of cluster models. Each family has the same initial concentration, from c = 0.8 for family a to 
c = 2.4 for family i. Letters correspond to the end models of each family. 
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Fig. 17. — The core collapse time as a function of the initial concentration for the models without tidal 
shocks. Filled dots are our results, triangles are from Quinlan (1996) , The agreement is very good, except 
for the very high concentration clusters. The solid line is our fit, eg. (|27|), with the last point excluded. 
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Fig. 18. — The core collapse time as a function of the shock parameter, p. The families of models with the 
initial concentrations c — 0.8 (a) and c = 1.0 (&) are excluded for clarity. Dotted lines are our fit, eq. (p8|). 
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Fig. 20. — The destruction time as a function of the shock parameter, /3. For clarity, some of the models are 
shown only for /3 < 0.03. Dotted lines are our fit, eq. (^. 



